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ABSTRACT 

To study how supernova feedback structures the turbulent interstellar medium, we construct 3D models 
of vertically stratified gas stirred by discrete supernova explosions, including vertical gravitational field and 
parametrized heating and cooling. The models reproduce many observed characteristics of the Galaxy such 
as global circulation of gas (i.e., galactic fountain) and the existence of cold dense clouds in the galactic disk. 
Global quantities of the model such as warm and hot gas filling factors in the midplane, mass fraction of ther- 
mally unstable gas, and the averaged vertical density profile are compared directly with existing observations, 
and shown to be broadly consistent. We find that energy injection occurs over a broad range of scales. There 
is no single effective driving scale, unlike the usual assumption for idealized models of incompressible turbu- 
lence. However, >90% of the total kinetic energy is contained in wavelengths shortward of 200 pc. The shape 
of the kinetic energy spectrum differs substantially from that of the velocity power spectrum, which implies 
that the velocity structure varies with the gas density. Velocity structure functions demonstrate that the phe- 
nomenological theory proposed by Boldyrev is applicable to the medium. We show that it can be misleading 
to predict physical properties such as the stellar initial mass function based on numerical simulations that do 
not include self-gravity of the gas. Even if all the gas in turbulently Jeans unstable regions in our simulation is 
assumed to collapse and form stars in local freefall times, the resulting total collapse rate is significantly lower 
than the value consistent with the input supernova rate. Supernova-driven turbulence inhibits star formation 
globally rather than triggering it. 

Subject headings: hydrodynamics — ISM: kinematics and dynamics — ISM: structure — methods: numerical 
— turbulence 



1. INTRODUCTION 

The interstellar medium (ISM) is observed to be turbulent. 
The electron density spectrum in the warm ionized medium 
shows a nearly Kolmogorov spectrum over at least six orders 
of magnitude in length (Armstrong et al. 1995). Supposing 
the small density perturbations measured to behave as a pas- 
sive scalar, this can be taken as evidence for a Kolmogorov 
velocity spectrum in that gas. Giant molecular clouds, which 
contain most of the molecular gas, show supra-thermal emis- 
sion linewidths (Zuckerman & Palmer 1974), indicating the 
presence of supersonic random motions. 

What drives this turbulence? Massive stars are a major 
structuring agent in the ISM (McCray & Snow 1979). Once 
formed inside molecular clouds, they influence the surround- 
ing medium via intense ionizing radiation and stellar winds 
until they end their lives in catastrophic explosions as super- 
novae, heating and stirring the ISM throughout their lives. 
Among these energetic processes, supernovae (and their col- 
lective effect, superbubbles) are likely to be the dominant 
contributor to the observed supersonic turbulence (Norman & 
Ferrara 1996; see Mac Low & Klessen 2004 for a recent re- 
view). This turbulence changes the balance between pressure 
and gravity on a wide range of scales, from galactic (Boulares 
& Cox 1990) to sub-cloud scales (e.g., Larson 1981), and in 
turn may regulate subsequent star formation, hence the term 
"feedback." With strong enough feedback, a galactic-scale 
superwind may be driven out of the galaxy (Heckman et al. 
2000,2001). 

Our understanding of how the radiative and hydrodynamic 
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feedback affects spatio-temporal structure of the ISM and star 
formation is only qualitative. It remains difficult to incorpo- 
rate the feedback into a consistent theory of star formation and 
the ISM, partly because star formation is an inherently com- 
plex phenomenon that involves an interplay between gravity, 
magnetohydrodynamics, and thermodynamics of compress- 
ible gas. Yet it is a crucial physical process for correct predic- 
tions of, e.g., the thickness, disk stability, and star formation 
rate of galaxies. Many of the serious problems that cosmo- 
logical models are faced with today are also believed to be 
associated with such feedback. 3 

Since no analytic theory based on first principles exists 
for turbulence in a compressible medium (see, however, a 
heuristic model by Sasao 1973), most previous work has re- 
lied heavily on numerical methods. Recent numerical sim- 
ulations in periodic boxes yielded key insights into the role 
of turbulence in star formation. First, supersonic turbulence 
in an isothermal or an adiabatic medium was found to de- 
cay quickly, i.e. within a sound crossing time, whether the 
medium was magnetized or not (Stone et al. 1998; Mac Low 
et al. 1998; Padoan & Nordlund 1999), unless it was con- 
stantly stirred by an external source of driving (Mac Low 
1999). Second, supersonic turbulence delays or sometimes 
prevents collapse (Klessen et al. 2000, hereafter KHM00), 
and is hence thought to be responsible for the observed low 
overall star formation efficiencies (Vazquez-Semadeni et al. 
2005). Overall, turbulence impedes collapse in unstable re- 
gions. Third, in contrast to the second point, it also creates 

3 The standard CDM paradigm predicts basic properties of the luminous 
component of matter on galactic and sub-galactic scales which are seriously 
inconsistent with observations. Outstanding astrophysical problems include 
incorrect predictions for the shape of the galaxy luminosity function and for 
the fraction of gas that should have cooled to form galaxies (the overcooling 
problem) (Somerville & Primack 1999; Cole et al. 2000). 



2 



Joung & Mac Low 



converging flows that may enhance density and thus promote 
collapse (Ballesteros-Paredes et al. 1999). Even when turbu- 
lence can support a cloud globally against gravitational col- 
lapse, it produces density enhancements that allow local col- 
lapse in both non-magnetized (KHMOO) and magnetized me- 
dia (Heitsch, Mac Low, & Klessen 2001, hereafter HMK01). 
However, the net effect of turbulence is to inhibit collapse 
(Mac Low & Klessen 2004). 

In some of these previous simulations, turbulence was 
seeded as the initial condition and left to decay for the rest 
of the run (Porter et al. 1992; Porter et al. 1998; Stone et al. 
1998; Mac Low et al. 1998; Padoan & Nordlund 1999). In 
other simulations, turbulence was driven constantly but the 
driving was included only in a general way as ubiquitous 
Fourier forcing, i.e., a random forcing in a prescribed nar- 
row range of wavenumbers applied uniformly to the box (Mac 
Low 1999; KHMOO; HMK01). These models only partly rep- 
resent the real ISM, where forced and decaying regimes may 
coexist (Avila-Reese & Vazquez-Semadeni 2001) and where 
an appropriate forcing function is probably broad-band (Nor- 
man & Ferrara 1996). It has been demonstrated that the slope 
of the power spectrum for the forcing function profoundly 
affects the statistics of turbulent fluctuations (Bonazzola et 
al. 1987; Vazquez-Semadeni & Gazol 1995; Biferale et al. 
2004). Hence, it is of great interest to determine the wave- 
length distribution of the kinetic energy in a more realistic 
medium driven by discrete physical space forcing. 

Our goal is to examine physical characteristics of the ISM 
driven by realistic, discrete physical space forcing in 3D. In 
particular, we study density and velocity structures and deter- 
mine the kinetic energy power spectrum. If, as some previous 
works suggest, the interaction of blast waves from supernovae 
structures the interstellar gas, we must resolve the individual 
explosions in order to correctly reproduce the main statistical 
properties of the ISM. For an accurate simulation of super- 
nova driving, a wide range of length scales from small clouds 
(~ 1 pc) to at least the "driving scale" (using the language de- 
rived from incompressible turbulence research) must be re- 
solved. In addition, to predict the cloud structure, the model 
must take into account appropriate heating and radiative cool- 
ing processes. 

There have been modeling efforts toward more realistic su- 
pernova driving by including its explosive nature. Rosen, 
Bregman, & Norman (1993), Rosen & Bregman (1995), and 
Wada & Norman (2001) performed two-dimensional hydro- 
dynamic simulations for the ISM in a galactic disk includ- 
ing feedback from supernovae and/or stellar winds. Results 
from two-dimensional models must be interpreted with cau- 
tion, since turbulence behaves differently in 2D than in 3D: 
in 2D, an inverse cascade of energy occurs towards large 
scales, while the enstrophy, the square of vorticity, cascades 
to small scales (Kraichnan 1967; Frisch 1995). More re- 
cently, three-dimensional models have been studied by several 
groups. These authors find that blast waves from supernovae 
sweep up the ISM into relatively thin shells that collide with 
one another to form cold dense clouds. However, none of the 
goals raised above has been adequately addressed. Korpi et 
al. (1999) used resolution of 10 pc, insufficient to follow the 
details of turbulent interactions, and only covered (0.5 kpc) 2 x 
2 kpc, although they did include magnetic field and galactic 
shear. Their main interest was to simulate a turbulent galactic 
dynamo. Slyz et al. (2005) included self-gravity of the gas but 
suffer from low spatial resolution (> 10 pc) and the absence of 
vertical density stratification. Kim et al. (2001) and Mac Low 



et al. (2005) also included magnetic field, but not galactic 
stratification. Avillez (2000) and Avillez & Berry (2001) used 
a stratified model similar to ours, but their studies focused on 
the properties of a Galactic fountain in the Milky Way. All 
the aforementioned 3D models employed minimum gas tem- 
peratures T m j„ = 100-310 K, so no cold gas could form. The 
difference between Avillez & Breitschwerdt (2004a, b) and 
the present work is more subtle, and will be discussed in § 
l3~2l 

We simulate a small patch of our Galaxy, (0.5 kpc) 2 in 
area. The computation box is elongated in z and extends 
from -5 kpc to +5 kpc to study the vertical structure. We 
choose the length scales in our simulations (2 pc < / < 10 
kpc) to lie between those of star-forming molecular clouds 
and those of large-scale galactic outflows (Heckman et al. 
2000, 2001). Our box size represents an optimal choice to 
improve our understanding of this enigmatic yet relatively un- 
explored regime. Expanding shells and their interactions are 
well-resolved at this scale. We attempt to build local, high- 
resolution models of the ISM based on which subgrid models 
for turbulent pressure can be developed, that will provide in- 
sight into how supernova feedback should be treated in global, 
cosmological simulations. We will continue to pursue this 
problem in a companion paper (M. Joung & M. Mac Low 
2005, in preparation; hereafter Paper II). Here, in Paper I, we 
describe our numerical model (§EJ and present the basic re- 
sults (§ [3}- In § |4] we explore where gravitational collapse 
would occur in the presence of explosion-driven turbulence. 
Finally, we summarize our findings and discuss their implica- 
tions in § |5] 

2. THE MODEL 
2.1. Basic Features 

Our simulations are performed using Flash (Fryxell et al. 
2000), an Eulerian astrophysical hydrodynamics code with 
adaptive mesh refinement (AMR) capability, developed by the 
Flash Center at the University of Chicago. It solves the Euler 
equations using the piecewise-parabolic method (Colella & 
Woodward 1984) to handle compressible flows with shocks. 
For parallelization, the Message-Passing Interface library is 
used; the AMR is handled by the PARAMESH library. 

We set up a stratified gas at 1.1 xlO 4 K initially in hy- 
drostatic equilibrium. The surface mass density of the gas, 
"Egas = 7.5 M© pc" 2 . The computation box contains a volume 
of (0.5 kpc) 2 x(10 kpc), elongated in the vertical direction. At 
1.95 pc resolution, it effectively contains 256 2 x5120 zones. 
We modify the refinement and derefinement criteria so that it 
becomes gradually more difficult to refine (and easier to dere- 
fine) with increasing distance from the midplane (z = 0). In 
practice, the inner ^600 pc near the midplane ends up refined 
to the maximum level. This allows us to capture the vertical 
flow, while not unnecessarily resolving low-density regions 
far from the disk. For our fiducial model, the number of dy- 
namically allocated zones is only 7-8 % of that in a single 
mesh code with the same maximum resolution. 

The gravitational potential of Kuijken & Gilmore (1989) is 
employed. It yields the vertical gravitational acceleration 

g(z) = a ' Z -a 2 z, (1) 



where a x = 1.42 X 10 -3 kpc Myr" 2 , a 2 = 5.49 x 10" 4 Myr" 2 , 
and zo = 0.18 kpc. The two terms on the RHS of the equation 
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represent contributions from a stellar disk and a spherical dark 
halo. It is static in time; self-gravity of the inhomogeneous, 
evolving gas is not included. Outflow boundary conditions are 
used on the upper and lower surfaces parallel to the Galactic 
plane, while periodic boundary conditions are used elsewhere. 

For each supernova explosion, we add thermal energy E m = 
10 51 ergs in a small sphere whose radius varies as a func- 
tion of the local density. The radii are chosen such that ra- 
diative losses inside the spheres are negligible in the first few 
timesteps after the explosions. The details are laid out in § 12.31 
below. We assume a constant overall rate of supernova explo- 
sions, and treat it as an input parameter. This assumption of 
constant supernova rate is justified by the fact that dynami- 
cally important variations in the supernova rate occur on ~1 
Gyr time scale, far longer than the duration of our simulations. 

In addition, each supernova produces ^500 metal particles 
that evolve passively following the gas component. We tag 
each particle with a unique number and follow its position, 
velocity, and thermal properties of the gas (density, tempera- 
ture, and pressure) in the cell that contains the particle. These 
variables are recorded at every timestep for all particles. This 
procedure enables us to trace the thermal history of metal par- 
ticles and to estimate their escape fraction from the galaxy. 

We ass ume that a fixed fraction (15/31 in the fiducial 
model; § 12. 3> of the supernovae are closely correlated in 
space as a way of simulating superbubbles. The remaining ex- 
plosions have random positions scattered through the galaxy 
mass to represent field supernovae. 

We solve the standard hydrodynamic equations, namely the 
continuity equation, the momentum equation, and the energy 
equation in a background gravitational field. Here we write 
only the energy equation in our model, which is modified to 
include supernova heating. In conservative form, 

^ + V-[(pE + P)v] = pvg-n 2 £ + S, (2) 
at 

where p, P, v, and T denote density, pressure, velocity, and 
temperature of the gas; the specific total energy E is the 
sum of internal energy e and kinetic energy per unit mass: 
E = e+ \v\ 2 /2. Gas pressure is related to density via the re- 
lation P = (7- l)pe with the adiabatic index 7 = 5/3. The 
gravitational acceleration g = g(z)z given by equation Q. The 
net cooling rate per unit volume n 2 C = n 2 A(T)-nT(z), where 
n is the number density of gas. For the cooling function, A(T), 
radiative cooling appropriate for an optically thin plasma with 
Z/Z© = 1 (cosmic abundance) is included assuming equilib- 
rium ionization (Sutherland & Dopita 1993). It is displayed 
in Figure [0 approximated as a piecewise power law. When 
the cooling function is expressed as A,(T) oc T 13 ', the gas 
is thermally stable where /3,- > 1. For the particular cool- 
ing curve that we adopted, the gas is thermally stable for 
10K<r<40Kandl0 4 K<r< 1.7 x 10 4 K, and marginally 
stable for 40 K < T < 200 K. The portion of the curve at low 
temperatures (T < 2 x 10 4 K) is adopted from Dalgarno & 
McCray (1972) assuming an ionization fraction of 10~ 2 . The 
cooling rate is most uncertain in this regime since the approx- 
imation of thermal equilibrium is usually invalid. Also shown 
in dotted line in Figure [^is the cooling curve from Spaans & 
Norman (1997), as reproduced in Wada & Norman (2001; see 
their Figure 1). The exact forms for the two heating terms in 
our model, the diffuse heating function T(z) and the impulsive 
he ating fr om lo cal supernova explosions S(x,t), are given in 
§§ !2.2l and l2.3l respectively. 
In dense regions, the cooling time f coo ; = (7- l)~ l kT /nA 
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FIG. 1. — Radiative cooling function for an optically thin plasma with 
Z/Zq = 1 approximated as a piecewise power law. The cooling curve is 
parametrized following Dalgarno & McCray (1972) with an ionization frac- 
tion of lfT 2 at T < 2 X 10 4 K and Sutherland & Dopita (1993) at T > 2 X 10 4 
K. Equilibrium ionization is assumed. For comparison, the cooling curve 
from Spaans & Norman (1997) is shown in dotted line. 

can be shorter than the Courant timestep Af. FLASH, as a 
default, applies heating and cooling terms in sequence us- 
ing an explicit method for updating temperatures. We have 
combined the separate heating and cooling routines so that 
temperature is updated only once in each timestep, and used 
an implicit method whenever the net cooling time t net = 
(7- \)~ x kT /nC < 0.1 Af. Thus, temperatures for high density 
(n > 10cm" 3 ) gas are computed using an iterative procedure, 
in our case, Brent's method (Press et al. 1992). 

2.2. Diffuse Heating 

It is crucial to choose the right diffuse heating rate. With- 
out this heating, high density gas will simply cool down to 
Tmin within a cooling time given by t coo \. In other words, no 
thermal equilibrium will exist above T m - m . Previous numeri- 
cal models adopted a diffuse heating rate T/,, which balances 
the vertical gravity g{z) on the gas to yield hydrostatic equi- 
librium at a uniform prescribed temperature 7}„„ (Mac Low 
et al. 1989; Avillez & Breitschwerdt 2004a). This heating 
rate is computed so that the net cooling rate in equation 10 is 
zero for all z. nT(z) = n 2 A(7/„,,), for the hydrostatic equilib- 
rium density profile given in equation[3]below. If the gas were 
isothermal, the cooling term would dominate where p > p/ 1s , 
while the diffuse heating term would dominate where p < p/ 1s . 

In contrast, photoelectric heating, i.e. photoemission of 
UV-irradiated dust grains (Bakes & Tielens 1994), has long 
been thought to be the dominant form of heating for the neu- 
tral component of the ISM, i.e., the cold neutral medium 
(CNM) and the warm neutral medium (WNM) (Wolfire et al. 
1995). It has the form T pe = 1.0 x lO-^eGoergs" 1 , where 
the heating efficiency e w 0.05 and Go is the incident far- 
ultraviolet field normalized to Habing's (1968) estimate of the 
local interstellar value (Gerritsen & Icke 1997). We adopt 
Draine (1978)'s value G = 1.7. 

Ionizing radiation from OB stars is an additional and a sig- 
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nificant heating source. Abbott (1982) estimates a total inte- 
grated luminosity e= 1.5 x 10~ 24 erg s" 1 cm" 3 in the disk of 
the Milky Way. Since only 20% of all available UV photons 
go into heating interstellar gas (Abbott 1982), the photoion- 
ization heating rate T p i w 3.5r pf . Much of this heating goes 
to the warm ionized medium. We do not, however, include it 
in our diffuse heating rate. 

Initially, the gas is assumed to be isothermal. We solve 
for the hydrostatic equilibrium profile of density for the given 
gravitational potential (see equation^: 

Phiz) = po exp I — — [^z 2 + zf ) -zo J - I , (3) 

where po and Po are the gas density and pressure at the mid- 
plane; Po = (Po/ fJfHn)kTo where 7b denotes the isothermal 
temperature of the gas. Note that the profile is narrower than 
a gaussian by the exponential factor due to a stellar disk (the 
term containing a{). This profile, ph s (z), would describe the 
vertical density distribution if the only form of support were 
thermal pressure. We evolve the systems for sufficiently long 
periods of time so that the memory of the initial condition is 
completely lost. 

Attempting to balance radiative cooling and diffuse heating 
of gas at our chosen initial gas temperature T^u = 1.1x10 K 
leads to I\, as 3QT pe . We believe that r/„, as used in previous 
models, is unrealistically high. It is the combination of proper 
(in fact, much lower) diffuse heating and supernova heating 
that together should balance radiative cooling. 

For the diffuse heating rate in our fiducial model, we adopt 
T pe , a value consistent with Wolfire et al. (1995). Note that 
our value of T is lower than that in Avillez & Breitschwerdt 
(2004a) by a factor of - 1 8 (their T m = 9 .0 x 1 3 K. The heat- 
ing rate coefficient T is assumed to be independent of gas den- 
sity; according to Wolfire et al. (1995), photoelectric heating 
weakly depends on density: T pe oc p 02 . 

The diffuse heating is applied to gas with T < T io „, where 
T wn = 2 x 10 4 K is the temperature at which hydrogen is fully 
thermally ionized. It is also assumed to decline exponentially 
in z as T oc ^oe~ ^ ' ^^,^ *. This exponential factor is a measure for 
the escaping ionizing radiation from the disk (Domgorgen & 
Mathis 1994; Dove & Shull 1994; Bland-Hawthorn & Mal- 
oney 1999, 2002), which should in principle be treated in a 
self-consistent way via a radiative transfer code (Fujita et al. 
2003; Wood et al. 2004). Existing observations are unable to 
provide a definite value of H^h for Milky Way type galaxies. 
In our model, Hdi, is set to 300 pc so that ^4 % of the ionizing 
radiation leaks out of the disk and into the halo, heating dif- 
fuse ionized gas there (Domgorgen & Mathis 1994). At high 
\z\, a minimum value r„„„ = 10~ 5 ro was used. Note that the 
prescribed functional form for T(z) precludes the possibility 
of obtaining thermal equilibrium for the initial gas distribu- 
tion at all heights. 

Our heating function V varies as a function of z only and 
is thus uniform at a given height in the disk (e.g., Gerritsen 
& Icke 1997). In reality, the diffuse heating rate is extremely 
non-uniform in space and time; it fluctuates by 1-2 orders of 
magnitude, and this large fluctuation in the heating rate may 
be responsible for the conversion of CNM to WNM and vice 
versa (Parravano et al. 2003; Wolfire et al. 2003). Thus, our 
use of constant diffuse heating rate is a simplification, justi- 
fied in part if non-uniformities created by far more energetic 
supernova explosions dominate the thermal energy budget of 
the gas. 



2.3. Adding Supernova Explosions 

One major obstacle that previous large-scale models includ- 
ing supernova feedback have faced was that when thermal en- 
ergy E sn = 10 51 erg was added locally into a sphere, most of 
the energy was radiated away too quickly. This cooling hap- 
pens within one or two timesteps after the explosion, i.e., be- 
fore blast waves form and start sweeping up the surrounding 
medium, converting the initially thermal energy into kinetic 
energy. As a result, adding impulsive heating in their models 
did not heat the medium in any significant way either ther- 
mally or dynamically (Katz 1992). To avoid this problem, 
each explosion should occur in such a way that no signifi- 
cant amount of energy is lost radiatively before the expand- 
ing sphere forms a strong blast wave and the interior density 
drops and relaxes to the Sedov profile. In order to achieve 
this effect, Thacker & Couchman (2001), for example, used 
an artificial time delay in cooling within the initial explosion 
sphere. Different algorithms such as assigning the explosion 
energy to fluid parcels as pure kinetic energy also have prob- 
lems (Navarro & White 1993). The fundamental cause of this 
problem is the lack of resolution in these models: relevant 
length scales for the energy source, in this case supernova 
explosions, lie below their resolution limits especially in the 
early phase of the expansion. In many previous models, an 
ad hoc thermalization efficiency of supernova energy, ranging 
from 10" 1 (Navarro & White 1993) to 10" 4 (Hernquist & Mi- 
hos 1995), was assumed. Because we track the thermal and 
dynamical evolution of blast waves from supernovae explic- 
itly, our model does not suffer from such limitations. 

Explosion radii in our model are determined so that the en- 
closed gas mass within each sphere M exp = 60M Q initially. In 
our fiducial model, this is less than 2 x 10~ 5 of the total gas 
mass in the computation box. We neglect the mass ejected by 
each supernova, i.e., no direct exchange of mass is assumed 
between stars and the interstellar gas. The radii usually vary 
between ^7 and ^50 pc. These initial spheres take up only 
a small fraction of the total volume. The mean local density 
inside the sphere is sufficiently low and the timestep is suf- 
ficiently short that no delay in cooling is necessary. Inside 
an explosion sphere, the density is redistributed uniformly to 
Pexp = 3M eA; ,/47rr 3 v , and then a thermal energy E m is injected 
evenly into the sphere. The timestep At w 10 3 yr. The value 
of Mexp represents a compromise in the sense that too small 
mass leads to too high cooling rates at T > 10 s K (due to 
bremsstrahlung) and too few zones within the explosion ra- 
dius. On the other hand, when the enclosed mass is too large, 
the initial cooling rates may be again too high if the injected 
thermal energy heats the sphere only to temperatures T < 10 6 
K. Also, the initial explosion sphere may occupy too much 
volume and thereby introduce artificial changes in the dynam- 
ics. 

As a preliminary test, we compared the time evolution of 
thermal and kinetic energies of a single supernova exploding 
in a uniform background medium (at various densities) with 
high resolution ID simulations by Cioffi et al. (1988). At 2 
pc resolution, their equation 3.15 provides a good fit for the 
energy evolution. 

The observed supernova frequency for the Galaxy is 1/330 
yr" 1 for Type I and 1/44 yr" 1 for Type II supernovae (Tam- 
mann et al. 1994). We normalize this to the area of our com- 
putation box so that the rates per unit area are 4.0 Myr" 1 kpc" 2 
and 30.0 Myr" 1 kpc" 2 . In the vertical direction, the frequency 
of supernovae is assumed to decline exponentially. The scale 
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heights for Type I and Type II SNe are 325 pc and 90 pc, re- 
spectively (Heiles 1987; Miller & Scalo 1979). 

It is known that spatial and temporal clustering of super- 
novae significantly affects their impact on the ISM (McCray 
& Snow 1979). However, most numerical simulations so far 
either completely ignored this correlation or treated clustered 
supernovae simply as one large constant luminosity wind bub- 
ble (Mac Low et al. 1989; Strickland & Stevens 2000). No- 
table exceptions are Korpi et al. (1999), who used statistical 
clustering of supernovae by introducing a density threshold 
for explosion sites, and Avillez (2000) and subsequent work, 
who assumed that 60% of all supernovae occurred at previ- 
ous explosion sites. No size distribution of superbubbles was 
specified by these models. Observations indicate that super- 
bubbles (or stellar clusters) follow a power-law distribution of 
the form dNg oc n~ 2 dn* where n» is the number of supernovae 
between a lower cutoff «„„■„ and an upper cutoff n max , and dN& 
is the number density of superbubbles having a total number 
of supernovae between «» and n*+dn* (Kennicutt et al. 1988; 
McKee & Williams 1997; Clarke & Oey 2002). Clarke & Oey 
(2002) point out that this distribution maximizes the mean 
volume occupied by each supernova. Interestingly, the n~ 2 
distribution places equal number of supernovae in any given 
decade in «». This is consistent with the large variations in 
the number of massive stars in clusters. 

For 3/5 of Type II supernovae, we explicitly account for 
a range of superbubble sizes in terms of the total number of 
supernovae that they contain, adopting the n~ 2 power-law dis- 
tribution. We fix «„„■„ = 7. To choose n max , we implement 
the following simple procedure. According to Kennicutt et 
al. (1988), there are 6-9 associations with 3500-7000 su- 
pernovae in the entire galaxy. Assuming the n~ 2 distribution, 
we choose n max so that there are about A ga i /At, ox associations 
with n 6 [n max /2, n max ] supernovae in the Galaxy. If we adopt 
Agai = tt(10 kpc) 2 and A\, ox = (0.5 kpc) 2 , n max w 40. In our 
model, all supernovae that belong to a particular superbubble 
explode at the same location, an approximation justified by 
the fact that once a stellar wind bubble forms, all subsequent 
supernovae explode inside it (Mac Low & McCray 1988). 
All superbubbles are assumed to have lifetimes of t SB = 40 
Myr, approximately the age of the least massive B star (8M Q ) 
when it explodes. We also assume that the given number of 
supernovae are equally spaced in time over t$a- During the 
first 5 Myr of each superbubble's lifetime, we include a con- 
stant mechanical luminosity wind originating from 4 pc ra- 
dius source region, mimicking stellar wind from the O stars. 
Shull & Saken (1995) argue that such early heating from stel- 
lar winds may accelerate shell growth and velocity evolution, 
compared to constant-luminosity models. The integrated lu- 
minosity of stellar wind heating is taken to be 0.14(n„£' s „) 
(Ferriere 1995). Mass outflows from stellar winds are not in- 
cluded. Isolated supernovae, which account for the remaining 
explosions, have random positions scattered throughout the 
disk. Since the progenitor stars drift apart with velocities ~5 
km s" 1 , they should lie out of their parent clouds by the time 
they explode, usually > 10 Myr (for Type lis) and > 1 Gyr 
(for Type la's) after they form. 

2.4. Missing Physics 

Aside from the limited spatial resolution of our simulations, 
our numerical model neglects several physical processes that 
are known to participate in interstellar cloud evolution and 
star formation. For example, since no self-gravity is included 
in our model, unstable regions do not collapse to form stars; 



instead, they just remain dense as a whole, so it is harder 
for shocks to destroy the clouds. A more complete model 
in the future should also include magnetic fields as well as 
shear from galactic rotation. Magnetic fields may substan- 
tially thicken swept-up shells (Ferriere et al. 1991) and pro- 
vide an additional pressure, which may the dominant com- 
ponent of the gas pressure at low temperatures (Heiles & 
Troland 2005; Avillez & Breitschwerdt 2004b). A large-scale 
shear due to differential rotation, when combined with weak 
magnetic fields, gives rise to the magnetorotational instability, 
which may dominate the turbulent velocity dispersion in the 
outer part of the disk (Sellwood & Balbus 1999; Dziourke- 
vitch et al. 2004; Piontek & Ostriker 2005). 

Another possibly important yet ignored physical process is 
thermal conduction. Koyama & Inutsuka (2004) argued for 
the importance of including thermal conduction by running 
two models with and without it, which showed discrepant 
global evolutions of quantities such as kinetic energy and to- 
tal number of dense clouds. They further argued that, to prop- 
erly include thermal conduction, the Field length Xf must be 
resoved. This is not feasible in our models, where the spa- 
tial resolution Sx^S> Xf- However, this crucial role of thermal 
conduction on cloud formation (via the action of thermal in- 
stability) may not carry over to the regime where the medium 
is externally driven. Based on 2D simulations that included 
stellar energy injection, Vazquez-Semadeni et al. (2000) sug- 
gested that stellar-like forcing erases the signature of the insta- 
bility, although their models also did not resolve Xf. Avillez 
& Mac Low (2002) and Avillez & Breitschwerdt (2004a) also 
argued that turbulent diffusion may be more important than 
thermal conduction for mixing gas of different temperatures. 

3. RESULTS 
3.1. Vertical Direction 

For the first few tens of megayears of the simulation, the 
heating from supernovae is insufficient to counter-balance the 
gravity. As a result, the disk initially cools and all gas col- 
lapses towards the midplane. After the collapse, shock waves 
bounce back and propagate away from the midplane in both 
directions, heating the gas at high altitudes. Almost identical 
behaviors were observed by Avillez (2000). As more super- 
novae explode, turbulent and thermal pressures build up in the 
disk. The total kinetic energy in the box exceeds the thermal 
energy at f w 15 Myr and thereafter. By t ps 30 Myr, enough 
supernovae have exploded that the blast waves from the mul- 
tiple explosions permeate the disk and their interactions reg- 
ulate the global gas structure. The average gas temperature 
rises steadily over this period. After t sa 60 Myr, the system 
reaches a statistical steady state, until the end of our simula- 
tion at f « 80 Myr. The entire run took ~8 x 10 4 CPU hours 
on the TeraScale Computing System at the Pittsburgh Super- 
computing Center. 

Figure |2ja, b, c) shows typical density, temperature, and 
pressure variations in the vertical direction, taken at t = 79.3 
Myr. A clumpy layer of cold dense clouds with a thickness 
of ^200 pc sits in the midplane. Above and below this layer, 
a galactic fountain (Shapiro & Field 1976) is set up, as ob- 
served in the Milky Way, enabling disk-halo interactions. In 
our model employing the Galactic supernova rate, the time- 
averaged mean mass flux through any horizontal surface ap- 
proaches zero, consistent with the negligible fraction of gas 
escaping the computation box (see below). 

The low density gas, which permeates the bulk of the vol- 
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FIG. 2. — Vertical slices of (a) density, (b) temperature, (c) pressure, (d) z component of the velocity, and (e) column density of gas (S,, = J pdy) at t = 79.3 
Myr. The large shells seen at z > 0.4 kpc are due to two high-altitude Type I supernovae that exploded 3.9 and 1.5 Myr ago. Note the presence of round cloudlets 
in (b) and (e) about 100-250 pc away from the Galactic plane. They are fragmented (super)shells and falling towards the midplane at several tens of km s . 
Filaments of neutral gas are found up to ~2 kpc away from the midplane. Note that only the inner 2 kpc of our 10 kpc vertical grid is shown. 



ume, moves away from the midplane with typical 



10 2 



km s . Some expanding superbubbles are launched with z- 
velocities reaching 300-400 km s" 1 at the bases. These large 
bubbles subsequently break out of the disk and directly dis- 
pose their thermal energy content to the halo. As the dense 
gas (predominantly fragments of supershells) rises, its tem- 
perature drops because of the decreasing diffuse heating rate 
as well as radiative cooling of the gas. As a result, neutral gas 
is present up to ^2 kpc above and below the midplane. 

Some of that gas later condenses and returns to the disk. 
This process occurs in several forms. First, cold (T w 10 2 
K in the core) and dense (n « 1 cm -3 ) cloudlets can be seen 
100-250 pc away from the midplane in Figure[2b, e). These 
cloudlets fall towards the midplane at velocities of a few to 
several tens of km s" 1 . They typically occupy a small vol- 
ume ^(10 pc) 3 , contain a mass 10-10 3 M Q , and are structured 
as a cold dense core plus a lower density tail, embedded in 
warmer gas. In addition, filaments of neutral clouds are found 
at larger heights (|z| w 0.25-2.0 kpc). With temperatures of 
order a few thousand degrees, and average densities n i=s 0.3 
cm -3 , several orders of magnitude higher than their surround- 
ings, these elongated filaments are estimated to contain 10 3 - 
10 M Q of neutral gas. Some extend beyond several hundred 
pc in length. Their sizes generally increase with height, but 
this may be an artificial effect caused by change in spatial res- 
olution. It is tempting to associate the objects within \z\ < 1 
kpc with intermediate velocity clouds (Wakker 2001). Neu- 
tral clouds at even larger heights (|z| > 1 kpc) may be iden- 
tified as the tangent point clouds found throughout the inner 
Galaxy (Lockman 2002). These clouds are thought to contain 




FIG. 3. — Vertical profiles of (a) gas density, (b) temperature, (c) pressure, 
and (d) metal density at t = 79.3 Myr averaged over x-y planes at constant 
heights. In (a), the dotted line represents our initial isothermal density profile, 
while the dot-dashed line shows the observed vertical profile of gas, which 
is the sum of molecular, neutral, and ionized gases. The inset displays an 
expanded view of the region near the midplane, \z\ < 1 kpc. 



a large fraction of the neutral gas in the Galactic halo. Given 
the crude approximations we make with regard to heating and 
cooling for the cold dense gas, we consider this a reasonably 
successful test of the model. 

At high altitudes (z > 1 kpc), extremely hot diffuse gas in 
the halo has very long cooling times t coo i(ri) = 51/n_3 Myr, 
where n_ 3 = «/(10 _3 cm- 3 ), if we take the cooling rate at 
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FIG. 4. — Cuts through the midplane (z = 0) showing distributions of (a) density, (b) temperature, and (c) pressure at ; = 79.3 Myr. Density and temperature both 
vary by about seven orders of magnitude. Because of an approximate inverse relationship between the two quantities, thermal pressure changes by 2-3 orders 
of magnitude. This much variation in pressure is inconsistent with ISM models based on full pressure equilibrium between phases. Correlated supernovae are 
responsible for producing most of the hot gas in the midplane. 



T = 10 6 K. Avillez & Breitschwerdt (2004a) ran a similar 
model for a much longer period (400 Myr) and showed that 
after t w 70 Myr the statistical properties related to the vertical 
distribution of gas do not change appreciably. 

Figure[3]shows how the gas density, temperature, pressure, 
and metal particle density averaged over the x-y plane vary in 
the vertical direction. For comparison, we show in the dot- 
dashed line the observed vertical profile of gas, which is the 
sum of three components: molecular (Clemens et al. 1988), 
neutral (Dickey & Lockman 1990), and ionized (Reynolds 
1991) gases. The average density near the galactic midplane 
in our model is higher than the observed value by 2-3, while 
the density is somewhat underpredicted at a few disk scale- 
heights (0.1 kpc < |z| < 0.5 kpc). Uncertainties in observa- 
tions are large; see, e.g., Figure 10 of Dickey & Lockman 
(1990). The discrepancy, if real, implies that supernova driv- 
ing alone cannot quite provide the necessary support in the 
vertical direction to explain the observed distribution of gas. 
Additional components of pressure, e.g. from the magnetic 
field and cosmic rays, are expected to contribute significantly 
(Boulares & Cox 1990). Compared to our initial isothermal 
density profile (dotted) in equation [3] a larger amount of gas 
is present at high altitudes (Jz| > 0.5 kpc) due to the presence 
of a galactic fountain (Fig. |2j. 

Since outflow boundary conditions are used at the top and 
bottom surfaces, once the gas escapes with v z > 0, it never re- 
turns to the box. For the Galactic supernova rate, only a negli- 
gible fraction of the gas (~0. 1 %) escapes the computation box 
after 79.3 Myr of evolution. None of the ^3x 10 5 metal par- 
ticles created during the simulation escapes the grid bound- 
aries at z = ±5 kpc. The vertical distribution of metal parti- 
cles peaks at the disk midplane, just as that of the gas density 
does. However, the relative concentration of the metal par- 
ticles is substantially less pronounced, i.e., an enhancement 
factor ~10 30 instead of ^10 43 for the gas density, implying 
that metals are mostly associated with the warm or the hot gas 



phases that have larger scale heights than the cold gas. 

For comparison, we ran a model with eight times the Galac- 
tic supernova rate and four times the gas column density, 
scaled so that the Kennicutt relation is approximately satis- 
fied for the computation box. Galactic outflows of mass and 
metals are observed in this model. We will discuss the vertical 
structure of the gas in models with higher supernova rates in 
Paper II. 

3.2. In the Midplane 

Figure|4]displays distributions of density, temperature, and 
pressure on a horizontal slice through the midplane (z = 0) at 
t = 79.3 Myr. As the standard blast wave theory (Ostriker & 
McKee 1988) predicts, when the age of a supernova remnant 
approaches the cooling time of the shock-heated gas, the outer 
part of the bubble undergoes thermal instability and devel- 
ops a dense spherical shell. These shells collide with one an- 
other to form clouds that are filamentary in shape. The clouds 
have characteristic lengths of several tens to hundred parsecs. 
Note that these cold dense clouds form even in the absence 
of self-gravity, as previous numerical simulations have found 
(Rosen et al. 1993; Rosen & Bregman 1995; Korpi et al. 
1999; Avillez 2000). They are formed due to the collec- 
tive effect of thermal instability and supersonic turbulence. 
These clouds have temperature lower than 50 K, H2 forma- 
tion times of —1 Myr (Hollenbach & McKee 1979), and esti- 
mated masses ranging from ~40 3 M Q to ~5x 10 5 M Q . Hence, 
they can be identified as molecular clouds or giant molecular 
complexes observed in the ISM. Although in reality many of 
these clouds contain small scale sub-structures, such details 
are unresolved in our model. An analysis of the velocity field 
reveals that turbulent flows within the clouds are trans-sonic 
or supersonic on the cloud scale (see Fig. HOt : Paper II). 

Comparing Figures @Ja) and (b) reveals an approximate 
inverse relationship between density and temperature. The 
highest gas density reaches ^10 3 cm" 3 , and the temperature 
ranges from 10K to ~10 7 K. Both of these quantities vary by 
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FIG. 5. — Phase diagram showing the distribution of thermal pressure and 
density in the galactic midplane computed from eight time slices spanning 7.0 
Myr of time interval (from 72.3 to 79.3 Myr). The grey scale shows the log- 
arithm of the occupied gas volume. A considerable amount of gas is found in 
the thermally unstable regime (200 < T < 1.0 X 10 4 K). The straight lines in 
the upper left quadrant of the diagram are caused by the adiabatic expansion 
(P oc p 1 ) of young supernova remnants. Hence, the lines have a slope equal 
to the adiabatic index of gas 7 = 5/3. Interiors of young supernova remnants 
are associated with low densities and high pressures (P/k > 10 4 cm" 3 K). 
Thermal pressures in dense clouds reach values about an order of magnitude 
higher than the mean pressure in the general ISM, even in the absence of 
self-gravity. 

more than six orders of magnitude in the midplane. If they had 
an exactly inverse relationship, the thermal pressure would 
have been constant. However, as shown in the pressure map 
(c), there are regions that clearly depart from pressure equilib- 
rium. High pressure regions are usually associated with inte- 
riors of newly formed supernovae, shock-compressed regions, 
or dense cloud cores, whereas low pressure regions are often 
interiors of old supernova remnants. Overall, thermal pressure 
varies by 2 to 3 orders of magnitude, confirming the results 
of previous works (Padoan et al. 1997b; Passot & Vazquez- 
Semadeni 1998; Mac Low et al. 2005). This much variation 
in pressure is inconsistent with ISM models based purely on 
pressure equilibrium between phases induced by thermal in- 
stability (Field et al. 1969). 

For a detailed look at the pressure distribution, we display 
in Figure|5]a scatter plot of pressure vs. density (a phase dia- 
gram). Two stable branches of the thermal equilibrium curve 
(T«10 4 K and 10K < T < 40 K) are densely populated. How- 
ever, other parts of the diagram are occupied by a dynamically 
determined continuum of densities and temperatures, rather 
than several discrete phases. This is because lower density 
gas simply does not have enough time to cool back down to 
its equilibrium temperature before it gets shock-heated again. 
If we assume that a supernova shocks, on average, an area 
^(100 pc) 2 in the midplane, for the Galactic supernova rate, 
the mean time between shock passages will be ^1 Myr. This 




1 1 1 2 1 3 1 4 1 5 1 6 1 7 1 o 8 
T (K) 

FIG. 6. — Fractions of mass contained in logarithmically spaced tempera- 
ture bins for gas within 500 pc of the midplane. The hatched region indicates 
the thermally unstable regime for our cooling curve. Note that the tempera- 
ture range is more extended than in Heiles (2001), where it ranged from 500 
to 5,000 K. Some fraction of the cold gas (X,„) having temperatures below 
~25 K is expected to be molecular hydrogen. 

is much longer than the cooling time, ~50 Myr, for the gas 
having n = 10~ 3 cm" 3 and T = 10 6 K. A new physical descrip- 
tion is required for this explosion-dominated medium (Mac 
Low & Klessen 2004). 

Remarkably, ^70% of the points lie in a narrow range of 
pressure (P/k) between 10 2 5 and 10 35 cm" 3 K. This implies 
that the low to intermediate density (n < 10 2 cm" 3 ) ISM on 
several hundred pc scales should be described by a nearly iso- 
baric equation of state, instead of an isothermal one, as typ- 
ically assumed in many recent works on galaxy formation. 
The near-isobaric character of the medium is reinforced by 
the distribution of turbulent pressure, as detailed in Paper II. 

Physical quantities in our model can be compared directly 
with existing observations. In particular, we focus on sev- 
eral global quantities measured via H I absorption line stud- 
ies (Heiles 2001; Heiles & Troland 2003). They found: (1) 
>48% of the WNM by mass lies in thermally unstable tem- 
perature range; (2) ^60% of all neutral hydrogen (again, by 
mass) is in the WNM; and (3) The warm gas occupies ~50% 
of the volume in the Galactic plane. 

The distribution of mass in terms of gas temperature is dis- 
played in Figure[6] Based on the shape of our cooling curve, 
T = 200 K separates the CNM from WNM. The WNM ranges 
from 200 K to 1 .7 x 10 4 K. Some fraction of the cold gas is ex- 
pected to be molecular hydrogen. According to Heiles (2001), 
most of the mass in the CNM is contained in gas with temper- 
atures between 25 K and 75 K. It is then reasonable to assume 
a large fraction of colder gas (T < 25 K) to be molecular hy- 
drogen. We denote this mass fraction by X m , and argue that 
X m is close to unity. Adopting these definitions, we find: (1) 
67% of the WNM by mass is in the thermally unstable regime; 
(2) 53% of all H I is contained in the WNM for X m =\ (43% 
for Z m =0.9); and (3) The warm gas occupies 43-52% (49% 
on average) of the volume in the midplane. Overall, the re- 
sults are broadly consistent with the observations despite our 
neglect of magnetic fields and self-gravity. 

Typical hot gas filling factors in the midplane of our model 
f h w 41-51% ((f h ) =44%; see Fig. |4p). These values of 
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FIG. 7. — Probability density functions for (a) gas density, (b) temperature, and (c) pressure of gas near the midplane (|z| < 125 pc), taken from eight snapshots 
spanning 7.0 Myr (from 72.3 to 79.3 Myr). The sizes of subboxes (kernels) used to construct the PDFs increase by factors of 2 from black (1.95 pc) to red (125 
pc). The temperature PDF (b) demonstrates the extent to which the classical three-phase medium description is valid. The pressure PDF (c) shows that there is 
as much volume below the average pressure as there is above it. 



fh are lower than the estimate by McKee & Ostriker (1977) 
(f h ps 70%). They ignored supernova clustering and the ver- 
tical density stratification, which are expected to reduce //,. 
For example, models that include superbubble evolution pre- 
dict ff, i=s 17-20% near the solar circle (Heiles 1990; Fer- 
riere 1998). The interstellar magnetic pressure will further 
decrease the volume occupied by hot gas by resisting the ex- 
pansion of bubbles (McKee 1990; Ferriere et al. 1991; Slavin 
& Cox 1993). Existing observations of the hot interstellar gas 
are too limited to constrain /},, but in external galaxies it is 
thought that f h < 0.5 (Dettmar 1992; Brinks & Bajaja 1986). 

Contrary to our result, Avillez & Breitschwerdt (2004a) 
obtained low //, « 0.2 at 1.25 pc resolution. Their model 
is discrepant from ours in at least the following three as- 
pects, each of which contributes to lower //,: (1) Their back- 
ground heating rate was higher by a factor of ~ 18; (2) Their 
supernova rate was lower by ~1.7; and (3) They allowed 
Type II supernovae to explode only in cold (T < 100K) and 
dense (n > 10cm" 3 ) regions, while we used random locations. 
Higher background heating rates effectively shift the thermal 
equilibrium curve (see Fig. [5} diagonally — upward and to the 
right. As a result, the stable branch of the curve at T < 40 K 
mostly lies outside the pressure range occupied by the bulk 
of the gas. Since the mean thermal pressure in the midplane 
does not change when T is increased, this leads to less gas 
in the cold phase and more gas in the warm phase. We ex- 
perimented by running two models with cooling curves deter- 
mined using ionization fraction of 10 _1 , as adopted by Avillez 
& Breitschwerdt (2004a). When we applied 4.5 times the dif- 
fuse heating rate of our fiducial model, //, w 50-60%. We 
then restarted the model at t = 68 Myr, and increased the dif- 
fuse heating rate to 18 times the fiducial value. After only 5 
Myr of evolution, the volume filling fraction of the hot gas 
dropped to 35%. Instead, the warm gas occupied 50% of the 
midplane area. However, we believe that this elevated level 
of diffuse heating is unphysical. Using density and tempera- 
ture thresholds for supernova locations lowers //, because (1) 
The supernovae then destroy massive clouds, returning the gas 
to the intercloud medium; and (2) Supernovae that occur in 
dense regions occupy smaller volumes when integrated over 
their lifetimes (Clarke & Oey 2002). The results in Korpi et 
al. (1999) can be interpreted using similar arguments. 

3.3. Density Fluctuations 



3.3.1. Probability Density Function 

There are numerous ways to characterize density distribu- 
tions of gas. The simplest and the easiest to interpret is the 
probability density function (PDF) for the gas density. It is 
displayed in Figure 0a). We use various subbox sizes over 
which density PDFs are computed. Purple corresponds to the 
smallest cubic subboxes with 3.91 pc on a side, and red to the 
largest subboxes with 125 pc on a side. 

We find that most of the simulation volume is occupied 
by low-density gas due to supernova feedback, in accord 
with Slyz et al. (2005), who simulated turbulent interstel- 
lar medium models in a non-stratified (1.28 kpc) 3 box with 
periodic boundary conditions and >10 pc resolution. Stellar 
feedback and/or self-gravity of gas were included in some of 
their models. In terms of the density PDF, they found that 
(1) The models without stellar feedback showed log-normal 
PDFs with a single peak; (2) A power-law tail developed in 
the high-density end when self-gravity was included; and (3) 
The PDF became markedly bimodal when stellar feedback 
was included, regardless of whether self-gravity was included 
or not. 

Although our PDF is not bimodal, there is a hint of a broad 
peak near n = 10 cm" 3 for the (thermally stable) cold high- 
density gas. It is controversial if this part of the PDF can 
be approximated by a log-normal function. Using 2D simu- 
lations including heating, cooling, and self-gravity, Wada & 
Norman (2001) claim the high-density end of their density 
PDFs is well-fitted by a log-normal function (see their Fig- 
ure 16). In contrast, several numerical experiments that in- 
cluded either self-gravity of gas or non-isothermal equation 
of state (7^1) reported that high-density tails develop in non- 
isothermal cases (Scalo et al. 1998; Li et al. 2003). Although 
a log-normal density PDF is a natural outcome for isothermal 
gas (Passot & Vazquez-Semadeni 1998), support for such a 
PDF for non-isothermal cases remains weak. 

The one-point density PDF does not contain information on 
how dense regions or voids are connected in space. For this 
reason, the cloud mass spectrum cannot be derived from the 
density PDF alone, as Scalo et al. (1998) point out. To sup- 
plement the density PDF, we attempt two methods of analysis. 
First, Figure 0a) shows how the density PDF changes as the 
subbox size increases. That the PDFs show markedly differ- 
ent shapes attests to the importance of specifying the smooth- 
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FIG. 8. — (a) Angle-averaged density power spectrum, displaying a wide 
peak around kL/2ir fa 20. The box size L = 0.5 kpc. The power falls off 
at large wavenumbers (small wavelengths) due to numerical diffusion, (b) 
Kinetic energy spectrum (solid) and angle-averaged velocity power spectrum 
(dotted). The kinetic energy is distributed over a wide range of wavenumbers. 
There is no single effective driving scale, unlike in Kolmogorov's idealized 
picture of incompressible turbulence. We find that 90% of the total kinetic 
energy is contained in wavelengths A < 190 pc. Because of the highly inter- 
mittent density structure, the velocity power spectrum is not parallel to the 
kinetic energy spectrum, especially at large scales. To guide the eye, two 
straight lines are plotted: the Kolmogorov energy spectrum (dashed) and the 
spectrum containing an equal amount of energy per decade (dot-dashed). 

ing scale with which the density PDF is measured. More sig- 
nificantly, the Jeans mass Mj and the density threshold for 
gravitational collapse /?,/, also change as a function of scale 
(see Fig. II 0i . Despite the fact that we do not fully under- 
stand what determines the shapes of the PDFs as the subbox 
size increases, this clearly has implications for gravitational 
collapse, as we explore in § |4] 

3.3.2. Power Spectrum 

The other method to characterize density fluctuations is to 
use second-order statistics such as the auto-correlation func- 
tion or its Fourier transform (FT), the power spectrum. Our 
computation box is periodic in the x and y directions but 
non-periodic in z, while FFTs assume periodicity in all di- 
rections. Hence, before taking FFTs, we apply to the z- 
components of the variables the Hanning window defined by 
w(z) = (l/2)(l + cos(27rz/L ; )) where L, is the vertical extent 
of the box used for the FFT (0.5 kpc) and -L z /2 <z< L z /2. 

While the velocity power spectrum in an incompressible 
medium |w(a;)| 2 has the clear physical meaning of specific ki- 
netic energy (see the discussion of Parseval's theorem in § 
I3.4> . p(x) 2 does not. If the density PDF turns out to be log- 
normal and if s = \ogp can be taken as a Gaussian random 
variable, a power spectrum of s can completely specify the 
density distribution. Even though this is not the case (see 
Fig- EtX the density power spectrum is nevertheless useful 
for comparison with previous numerical models. 

The density power spectrum of an explosion-driven, 
strongly compressible medium is displayed in Figure[SIa). Its 
shape contrasts drastically with its counterpart in a weakly 
compressible medium, where density fluctuations behave as 
a tracer field and possess a Kolmogorov spectrum (Lithwick 
& Goldreich 2001). The spectrum peaks near the smallest 
scale at which kinetic motions can be well-resolved A p « 20 
pc. This wavelength is slightly smaller than the most energy- 
containing scales Xe ~ 20-40 pc (Fig.[8j>). Since X p is within 
the range affected by numerical diffusion (§ I3.4> . this conclu- 
sion is somewhat uncertain. We do note that our value of A p is 
comparable to the correlation length of 28 pc measured from 
two-dimensional autocorrelation functions of the molecular 



hydrogen column density map in the Taurus molecular com- 
plex (Kleiner & Dickman 1984). 

It is simple to predict the density power spectrum of a 
medium dominated by shocks, because a shock is essentially 
a step function in density and velocity fields. Step functions 
are generally associated with power spectra having P(k) ex k~ 2 
(Passot et al. 1988). However, neither observations nor nu- 
merical simulations support this relation for density. Padoan 
et al. (1997a) showed that the density power spectrum has 
a power-law form P(k) oc k~ a with a = 2.6 ± 0.5, similar to 
values determined from observations (Stutzki et al. 1998). 
However, comparison with observations is not straightfor- 
ward, since only the projected (column) density, not the phys- 
ical density, is available from observations. The steep spec- 
tral slope at larger wave num bers may have been caused by 
the nature of PPM; see § 13.41 Unlike the spectrum in Padoan 
et al. (1997a) who used a periodic isothermal box, our den- 
sity spectrum turns over at kL/ln ss 10 and falls off at small 
wavenumbers (large wavelengths). 

Finally, we add a cautionary note. Generally, one should be 
careful when computing power spectra of physical quantities 
from an AMR code. This is because power at small scales 
can be easily underestimated in relatively under-refined re- 
gions. To avoid such mistake, the choice of refinement and 
derefinement criteria is crucial (Kritsuk et al. 2004). Al- 
though our simulations are performed with an AMR code, the 
region near the midplane, particularly the cube with volume 
(500 pc) 3 within |z| < 250 pc, ends up refined uniformly to 
the maximum refinement level after a few tens of Myr of evo- 
lution due to extreme density inhomogeneities and turbulent 
motions caused by the explosions. (For this reason, in fact, we 
effectively turn off the AMR machinery at late times.) Thus, 
at late times, we may compute power spectra without misrep- 
resenting them. 

3.4. Energy Fluctuations: Driving Scale 

How is the kinetic energy distributed among various length 
scales in our turbulent medium stirred by supernova explo- 
sions? This question is crucial for star formation. Depending 
on the slope of the forcing function k~ a , the statistics of turbu- 
lent fluctuations and clumps change (Bonazzola et al. 1987; 
Biferale et al. 2004). 

Previous simulations of uniformly driven (as opposed to 
explosion-driven) turbulence with variable driving wave- 
lengths have shown that the occurrence and efficiency of local 
collapse into dense cores — and, presumably, stars — decreases 
as the driving wavelength decreases (KHM00). Despite its 
relevance to cloud formation, as far as we know, the driving 
scale has never been measured in numerical simulations of 
ISM driven by point explosions. Observations of molecular 
clouds show that power-law scaling extends up to the largest 
observed clouds (Ossenkopf & Mac Low 2002), and suggest 
that the turbulence is driven by large-scale, external sources. 
However, gas density and velocity — two physical quantities 
relevant to the kinetic energy density — can only be disentan- 
gled by observations under certain conditions (Brunt & Mac 
Low 2004; Lazarian & Esquivel 2003). 

In an incompressible or weakly compressible medium, 
gas density is nearly uniform and thus the energy spec- 
trum is simply the angle average of the velocity power spec- 
trum: E inc (k) = 4nk 2 (P v (k))n where k = \k\ and P r (k) is the 
three-dimensional Fourier transform of each component of 
v, squared and then summed. The velocity power spec- 
tra of compressible interstellar gas have been studied in 2D 
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(Dahlburg et al. 1990; Vazquez-Semadeni et al. 1996; Wada 
& Norman 2001) and in 3D (e.g., Kritsuk et al. 2004). How- 
ever, none of the previous simulations performed in 3D dealt 
explicitly with power spectra of the medium driven by discrete 
physical space forcing. 

In a strongly compressible medium with high Mach num- 
bers, the gas density p is highly intermittent and can vary by 
several orders of magnitude. To measure how much kinetic 
energy is contained near a given wavenumber k, we should 
compute the power spectrum of ^fpv instead of v, since ac- 
cording to Parseval's theorem: 

r L/2 



■L/2 



[f( X )] d X : 



\f(k)\ 



dk . 



(4) 



where f(x) and f(k) are Fourier transform pairs, and L is the 
size of the periodic box. If we substitute V((x) for f(x) (where 
i = x, y, z), we obtain a relation for the incompressible case, 
which justifies the use of the velocity power spectrum for the 
kinetic energy spectrum. If we now substitute -y/pv, in place 
of f(x), we obtain the corresponding relation for a compress- 
ible medium, which relates the kinetic energy density (pvf) 
to the Fourier transform of y/pv,. As we demonstrate below, 
the power spectrum of ^fpv shows a markedly different shape 
from that of v. This point was briefly noted in the discussion 
of Padoan & Nordlund (1999). 

Figure|SIb) shows the one-dimensional kinetic energy spec- 
trum E^ n (k) = Ank 2 (P w (k)} $i where w = ^fpv. The measured 
slope is far shallower than the Kolmogorov spectral index 
of -5/3 for incompressible turbulence. In a strongly com- 
pressible, explosion-driven turbulence, shocks and interact- 
ing blast waves contribute to a wide range of driving wave- 
lengths. The shallow spectral slope shows the presence of 
substantial power at small length scales, suggesting that SN- 
driven motions dominate the kinetic energy spectrum even at 
small scales where H II regions (Haverkorn et al. 2004) or 
protostellar outflows may also be active. Thus, turbulence in 
the dense medium is at least partly driven by supernovae. 

The angle-averaged velocity power spectrum often used in 
the literature is also plotted in Figure |Ub) for comparison. 
This spectrum contains velocity information only and neglects 
the density inhomogeneity, as though the medium were in- 
compressible. The remarkable difference at low wavenum- 
bers between the two spectra in Figure |8jb) suggests that 
dense regions, which more or less determine the shape of the 
kinetic energy spectrum, have a vastly different velocity struc- 
ture from rarefied (intercloud) regions. 

Although kinetic energy is distributed over a broad range of 
scales, 90% of the total kinetic energy is contained on wave- 
lengths A < 190 pc. This value (190 pc) is likely to decrease 
somewhat in higher resolution simulations that better resolve 
turbulent motions on < 20 pc scales. In typical cosmologi- 
cal simulations, therefore, most of the turbulent energy due to 
multiple supernova explosions is contained in scales smaller 
than a resolution element. To include the effect of hydrody- 
namic feedback in such simulations, it is necessary to use 
a subgrid model that represents unresolved, small-scale mo- 
tions (Paper II). 

To test the robustness of our results on numerical resolution, 
we ran our fiducial model with the same box size at three dif- 
ferent resolutions (8, 4, and 2 pc). The energy power spectra 
for the three runs display identical shapes in the dissipation 
range, but are shifted horizontally with respect to one another 
by ~2. This occurs because, in a compressible medium, struc- 
tures on scales smaller than ^10 grid zones are affected by 



numerical dissipation (Porter et al. 1992; Mac Low & Os- 
senkopf 2000), determining the inertial range of the turbu- 
lent cascade. Finite difference algorithms directly damp ki- 
netic motions on small scales and convert their kinetic energy 
into heat. It follows that velocity differences on small scales 
(< lOAx) are somewhat underestimated. However, that may 
not affect the overall dynamics since thermal pressure should 
dominate the turbulent pressure on such small scales. 

The power spectra tend to flatten slightly near the 
wavenumber of maximum dissipation. This distinctive shape 
of the spectrum in the dissipation range (< 10 Ax) is caused 
by the bottleneck effect known to occur in hydro codes using 
PPM (Falkovich 1994; Porter et al. 1992; Porter et al. 1998). 

3.5. Velocity Structure Function 

The velocity power spectrum is steeper than the slope ex- 
pected for Kolmogorov turbulence (-5/3). In addition, since 
the spectral slope varies depending on the wavenumber, it is 
unclear whether an inertial range exists at all. How can we 
understand this behavior? Extensive experimental and nu- 
merical studies have found deviations from the Kolmogorov 
spectrum, usually referred to as intermittency corrections. To 
quantify deviations from the Kolmogorov spectrum, it is use- 
ful to compute a velocity structure function defined by 



S p (l)= (\v(x + l)-v(x)\ p ) <xl C(p) , 



(5) 



where the scalar lag I = \l\. There are several phenomenolog- 
ical frameworks for compressible turbulence, each of which 
gives a distinct prediction for the scaling exponents £(p). We 
apply the idea of extended self-similarity, which showed cor- 
rect scaling behaviors for £(p) even in systems with moderate 
Reynolds numbers, specifically the power index -(1 + C(2)) = 
-1.74 for a relatively wide range in Si(l) (Benzi et al. 1993; 
Camussi & Benzi 1996). We plot the structure function 
against the third-order structure function Si(l) in Figure|9ja). 

Predictions for ((p) from various analytic theories of tur- 
bulence are shown in Figure |9jb). Based on She & Lev- 
eque (1994)'s model of intermittency for incompressible tur- 
bulence (see also She & Waymire 1995), Boldyrev (2002) 
proposed an extension, noting that the most disspative struc- 
tures in strongly compressible ISM tend to be shocks, i.e. 
two-dimensional sheets. Our structure functions for the trans- 
verse component of the velocity field are most consistent with 
the prediction from this theory. Those for the longitudinal 
component diverge from straight lines at large S3U) (i.e., large 
separations), because most of the gas in terms of volume trav- 
els away from the midplane. If we omit the last two data 
points for each value of p, we find that the longitudinal struc- 
ture functions also agree with Boldyrev (2002)'s prediction. 
Until our work, the theory has been compared well only with 
isothermal gas driven on large scales by a random Fourier 
space forcing (Boldyrev et al. 2002; Padoan et al. 2003). 
Our result demonstrates that the theory is applicable even to 
a medium driven by discrete point explosions and subject to 
nonlinear heating and cooling processes. 

4. GRAVITATIONALLY UNSTABLE CLOUDS 

The connection between small-scale star formation in local 
regions of cold dense gas and large-scale, galactic star forma- 
tion is unclear. The Kennicutt-Schmidt law is successful in 
predicting star-forming rates (SFRs) of galaxies with wildly 
different characteristics — from low-surface brightness galax- 
ies to starbursts — implying that the SFR can be estimated 
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FIG. 9. — (a) Longitudinal and transverse velocity structure functions, (b) 
Scaling exponents C(/>)- Theoretically expected scalings are shown for com- 
parison. From top to bottom, they are from Kolmogorov (1941) (dotted), She 
& Leveque (1994) (dot-dashed), and Boldyrev (2002) (dashed). The dot-dot- 
dashed line is plotted assuming D = 2.3, taken from Elmegreen & Falgarone 
(1996)'s observation of molecular clouds. The data points are from, again 
from top to bottom: longitudinal and transverse velocity structure functions. 



from macroscopic properties of galaxies, and is perhaps con- 
trolled mainly by gravitational instability. Such a claim was 
made by Li et al. (2005), who used an isothermal equation 
of state for gas with a high speed of sound to represent the 
interstellar turbulent pressure. One might argue, though, that 
star formation is inherently a local phenomenon that depends 
sensitively on the physical conditions on cloud scales. 

On kiloparsec scales, much larger than the driv ing s cales of 
the interstellar turbulence that we measured in § 13.41 gravity 
is important and the effect of turbulence can be represented 
as an effective turbulent pressure term (Elmegreen & Scalo 
2004). Below those scales, at the scales resolved by our nu- 
merical model, supersonic turbulence controls much of the 
cloudy structure. 

Observational and numerical studies have established that 
the Toomre (1964) parameter Q = na / '7rG£ gflJ , where k is 
the epicyclic frequency and a is the effective sound speed, 
controls whether or not a given region of a galaxy will ac- 
tively form stars (Klessen 1998; Li et al. 2005). However, 
even in Toomre-unstable regions, not all the gas collapses 
within the local freefall time. Several authors have argued that 
only a fixed (constant) fraction of mass fu above some den- 
sity threshold will collapse to form stars (Elmegreen 2002; 
Krumholz & McKee 2005), based on the universality of den- 
sity PDFs in turbulent media (Wada & Norman 2001). Still, 
many questions remain unanswered: what determines /m\ is 
there a clear density threshold for star formation; and is /m 
the same (on average) for various galaxy types? 

One of the aims of our work has been to study how gas 
dynamics and thermodynamics shape the ISM. In our model, 
cold dense clouds form directly in the turbulent flow. What 
fraction of these clouds are gravitationally unstable, given the 



density and velocity structures in our model? If gravity is 
decoupled from other cloud-forming processes such as tur- 
bulence and thermal instability, do we still reproduce a star- 
formation rate consistent with the value implicitly set by the 
input supernova rate? To answer these questions, we apply 
to the data a simple criterion for gravitational collapse (see 
below). By comparing the collapse rate with our input su- 
pernova rate, we can quantitatively check the validity of our 
high-resolution simulations as a model for star formation. 

Since our model does not include self-gravity of gas, we 
estimate the gravitational collapse rate (hence the SFR) by 
using a crude method that involves the modified Jeans crite- 
rion (Chandrasekhar 1951). This approach will yield the cor- 
rect result only if (1) the spatial resolution of our model is 
sufficiently high that the gas dynamics is well resolved, par- 
ticularly physical quantities such as mass density p and veloc- 
ity dispersion of gas a that are crucial for estimating the col- 
lapse rate, and (2) supersonic random motions determine most 
properties of star-forming molecular clouds, with just minor 
modifications from gravity (Padoan et al. 1997a; Padoan et al. 
1997b; Padoan et al. 1998). If self-gravity extensively alters 
the structure of clouds presented in § [3] our estimate of the 
collapse rate will be inaccurate. 

We analyze the gas in the rectangular region near the mid- 
plane within |z| < 125 pc that contains the bulk (>90%) of 
the total gas mass. We estimate the total collapse rate of gas 
as follows. First, we identify Jeans-unstable boxes at vary- 
ing scales using the criterion Mt ox /Mj > 1, where Mj = pXj, 
p is the average density in the box, Ay = (n /Gp) l l 2 <7 t ot, and 
a m = (c 2 + \a 2 ) 1 ! 2 (Chandrasekhar 1951). We use c 2 = 
J box c 2 dM and a 2 = (l/M box ) J box <7 2 dM. (Paper II 
contains more details on how c s and a are computed.) The to- 
tal velocity dispersion <j tot changes with scale, bec ause kinetic 
energy is distributed over a wide range of scales (§ 13. 4> . Then 
the total collapse rate is estimated assuming that the unstable 
boxes collapse within the local freefall time. To convert this 
value to the SFR, we multiply it by a star formation efficiency 
of 0.30, as it is thought that ^30% of the cloud mass that grav- 
itationally collapses is converted to stars (Li et al. 2005 and 
references therein). 

Figure fTol' a) shows that there is no single density threshold 
Pth for gravitational collapse. Instead, p,h depends on the scale 
at which collapse occurs. The turbulent velocity dispersions 
for gas with number density of 10 2 cm" 3 range from 0. 1 to 1 .0 
km s" 1 (Fig. 110b). For the small (purple and blue) subboxes, 
the turbulent velocity dispersion is lower than the local sound 
speed (Fig. 1 10b). as turbulent motions on those scales are not 
well resolved in our model. 

The plot of the predicted SFR against the subbox sizes is 
displayed in Figure Our results show that, for L box > 20 
pc where turbulent motions are resolved, even if all the gas 
in turbulent Jeans unstable boxes collapses and forms stars 
within the local freefall time f/y, the resulting star formation 
rate remains lower than the value consistent with our input su- 
pernova rate: 4.4 (or 6.8) x 10~ 3 M yr _1 kpc~ 2 , which is com- 
puted assuming 130 (or 200) M Q of stellar mass is required 
per supernova using the Salpeter IMF. 

Turbulent compression and the structuring of the ISM by 
supernovae may set the stage for further gravitational collapse 
and even determine where it occurs. However, the discrep- 
ancy between the assumed and the predicted SFRs illustrates 
a limitation of our model. We interpret this as an indication 
that the density PDF must be affected by self-gravity; this idea 
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FIG. 10. — (a) Mass contained in each subbox vs. average gas density of the subbox. There is no single density threshold for gravitational collapse, (b) Turbulent 
velocity dispersion vs. average density of gas in the subbox. (c) Same as (b), but plotted against the mass-weighted gas temperature. Most Jeans-unstable boxes 
on the cloud scale are associated with transonic or supersonic velocity dispersions. 




10 100 
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FIG. 1 1 . — Predicted star formation rate from the model plotted against the 
subbox sizes used. The dotted line is drawn assuming that 30% of the mass in 
Jeans unstable regions turns into stars. The red dashed lines show the SFRs 
consistent with the assumed Galactic supernova rate, assuming 130 or 200 
Mq of stellar mass is required per supernova. 

is supported by numerical results that the high-density end of 
the PDF, most crucial for gravitational collapse, is positively 
skewed when self-gravity is turned on (Li et al. 2003; Slyz 
et al. 2005). If this is true, computations of SFRs based di- 
rectly on the assumption of a log-normal density PDF may be 
misleading. Uncertainties in the present analysis are missing 
kinetic energy caused by the finite resolution (§ I3.4> . which 
tends to increase predicted SFRs on small scales, and the ab- 
sence of magnetic fields. Magnetic pressure reduces the prob- 
ability of collapse (HMK01; Vazquez-Semadeni et al. 2005). 

We conclude that, while it may be possible to predict where 
gravitationally bound regions will occur, estimating the SFR 
based on our model leads to an underestimate, about an or- 
der of magnitude lower than the value consistent with the 
assumed supernova rate. This limitation of our ISM model 
probably stems from the key physical process that we ignore: 
self-gravity of the gas. Other models relying entirely on tur- 
bulent compression will likely show similar deficits. 

We verified that the gas sitting at a uniform temperature 
(T = 4.2 x 10 3 K, corresponding to the total (thermal plus tur- 
bulent) velocity dispersion of our medium a ~ 6.5 km s" 1 ) in 
hydrostatic equilibrium with the imposed gravitational field is 
stable. In other words, no Jeans unstable boxes occur. How- 
ever, this point cannot be used to argue that supernovae en- 



hance the global gravitational collapse rate and trigger star 
formation. In fact, the opposite is true. Since the gas was ini- 
tially out of thermal equilibrium, it would have promptly col- 
lapsed and cooled into a thin sheet. Most of the mass would 
have become gravitationally unstable, without extra stirring 
from supernovae. The net effect of supernova-driven turbu- 
lence is to inhibit star formation globally by decreasing the 
amount of mass unstable to gravitational collapse. 

5. SUMMARY 

In this work, we have constructed 3D models of the ISM in- 
cluding vertical stratification and discrete physical space forc- 
ing by supernova explosions. We included parametrized heat- 
ing and cooling as well as both isolated and correlated super- 
novae. This work supports models that treat supernova-driven 
turbulence as a major structuring agent in the ISM. Our main 
results are as follows: 

• Cold dense clouds naturally form as blast waves sweep 
the gas and collide. The intercloud region is filled 
with warm or hot gas having relatively low densities. 
The resulting density and velocity structures adequately 
match the observed filamentary and clumpy structure of 
the ISM. 

• A galactic fountain reaching a height of several kilo- 
parsecs forms, with most of the volume moving away 
from the midplane. Fragmented shells rain back down 
as small cold clumps, which may be identified as inter- 
mediate velocity clouds. 

• Global quantities in our model are broadly consistent 
with existing observations. Specifically, our results are 
consistent with observations of the mass fraction of 
thermally unstable gas and the warm gas filling factor 
in the Galactic plane derived from H I absorption lines 
(Heiles 2001 ; Heiles & Troland 2003). 

• The diffuse heating rate T influences the global struc- 
ture of our model. For example, in our numerical exper- 
iments, increasing T by 4 lowered the occupation frac- 
tion of the hot gas in the Galactic plane by ~20%. This 
occurs because high background heating rates effec- 
tively shift the thermal equilibrium curve in the phase 
diagram so that the cold (thermally stable) branch lies 
mostly outside the range of pressure occupied by the 
bulk of the gas. 
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• The vertical distribution of gas in our model deviates 
from the profile inferred from observations. Specifi- 
cally, the gas is excessively concentrated in the disk 
midplane by factors of 2-3, while it is deficient at 
01 ^5 \z\ < 0.5 kpc. We attribute this discrepancy to 
neglected components of pressure such as those from 
the magnetic field and cosmic rays. 

• The density power spectrum has a broad peak close to 
the dissipation scale A ~ 20 pc. It turns over and falls 
off at large scales. 

• There is no single effective driving scale; instead, we 
see energy injection over a range of scales. As a result, 
the kinetic energy power spectrum possesses a slope 
much flatter than Kolmogorov's spectrum at scales be- 
tween ^20 and ^500 pc. The predominance of small 
scale density structure causes the kinetic energy power 
spectrum to look remarkably different from the often- 
used velocity power spectrum. About 90% of the total 
kinetic energy is contained in wavelengths shortward of 
190 pc. 

• The velocity structure functions demonstrate that the 
phenomenological theory proposed by Boldyrev (2002) 
is applicable even to a medium driven by discrete point 
explosions and subject to nonlinear heating and cooling 
processes. 

• Even if all the gas in turbulent Jeans unstable subboxes 
in our simulation collapses and forms stars in local 
freefall times, the resulting collapse rate is significantly 
lower than the value consistent with the input super- 
nova rate. Supernova-driven turbulence inhibits, but 
does not prevent, star formation. From this, we infer 
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that the density PDF must be significantly affected by 
self-gravity. For predicting physical properties of star- 
forming regions (such as the stellar initial mass func- 
tion), the statistics of turbulent fluctuations generated 
by models without self-gravity can be misleading. 



Applying high supernova rates to our model galaxy nat- 
urally leads to galactic outflows that have already been ob- 
served at both low (Heckman et al. 2000) and high redshifts 
(Pettini et al. 2001, 2002; Shapley et al. 2003). In the future, 
we will use these local models to measure the efficiency of 
energy transfer into superwinds out of a stratified disk for a 
given star formation rate, and thereby construct subgrid mod- 
els for stellar energy feedback in cosmological simulations. 
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